{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Factorizaciones y otras diversiones\n",
    "Autor: Andreas Noack Jensen (MIT) (http://www.econ.ku.dk/phdstudent/noack/)\n",
    "(con edición de Jane Herriman=\n",
    "\n",
    "## Sinopsis\n",
    " - Factorizaciones\n",
    " - Estructuras de matrices especiales\n",
    " - Álgebra lineal genérica"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Empezamos con un sistema lineal de la forma\n",
    "\n",
    "`Ax = b`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Array{Float64,2}:\n",
       "  0.157959  0.838031  -0.434173\n",
       " -0.404891  0.306374   0.79616 \n",
       "  0.77716   2.11739   -0.28197 "
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = randn(3,3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Array{Float64,1}:\n",
       " 0.561817\n",
       " 0.697643\n",
       " 2.61258 "
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = fill(1.0, (3))\n",
    "b = A*x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Factorización\n",
    "La función `\\` esconde cómo regularmente se resuelve el problema.\n",
    "\n",
    "Dependiendo de las dimensiones de `A`, distintos métodos son elegidos para resolver el problema.\n",
    "\n",
    "```\n",
    "Ax = b\n",
    "```\n",
    "\n",
    "Un paso intermedio en la solución es el cálculo de la factorización de la matriz `A`. \n",
    "\n",
    "Básicamente, una factorización de `A` es una manera de expresar `A` como el producto de matrices triangulares, unitarias, y de permutación.\n",
    "\n",
    "Julia guarda estas factorizaciones usando un tipo abstracto llamado `Factorization` y varios subtipos.\n",
    "\n",
    "Un objeto `Factorization` entonces debería ser pensado como una represenatación de la matriz `A`.\n",
    "\n",
    "#### LU"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cuando `A` es cuadrada, el sistema linear es resuelto factorizando `A` vía\n",
    "\n",
    "```\n",
    "PA=LU\n",
    "``` \n",
    "\n",
    "donde `P` es una matriz de permutación, `L` es un triangular inferior unitaria y `U` es triangular superior. \n",
    "\n",
    "Julia permita calcular la facorización LU y define un tipo de factorización compuesta para guardarlo."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Podemos hacer una factorización LU sobre `A` vía `lu(A)` ó `lufact(A)`.\n",
    "\n",
    "La función `lu` regresa matrices `l` y `u` y un vector de permutación `p`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "([1.0 0.0 0.0; -0.520989 1.0 0.0; 0.203251 0.289228 1.0], [0.77716 2.11739 -0.28197; 0.0 1.40951 0.649257; 0.0 0.0 -0.564645], [3, 2, 1])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "l,u,p = lu(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "El pivoteo está prendido por default, o sea que no podemos assumir que A == LU.\n",
    "Vamos a comprobar esto viendo la norma de `LU - A`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.021547018129227"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "norm(l*u - A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Esto nos muestra que queremos tomar en cuenta el pivoteo!\n",
    "\n",
    "Podemos pensar en `A[p,:]` como la sintaxis para `PA`, o como el producto de una matriz de permutación con `A`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.6653345369377348e-16"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "norm(l*u - A[p,:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Por otro lado, podemos apagar el pivoteo usando `Val{false}` en Julia 0.6 ó `Val(false)` en versiones más modernas."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "([1.0 0.0 0.0; -2.56327 1.0 0.0; 4.92001 -0.817175 1.0], [0.157959 0.838031 -0.434173; 0.0 2.45447 -0.316742; 0.0 0.0 1.59533], [1, 2, 3])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "l,u,p = lu(A, Val{false})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cuando apagamos el pivoteo, `LU = A`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.1102230246251565e-16"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "norm(l*u - A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Una segunda manera de hacer la factorización es con `lufact`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Base.LinAlg.LU{Float64,Array{Float64,2}} with factors L and U:\n",
       "[1.0 0.0 0.0; -0.520989 1.0 0.0; 0.203251 0.289228 1.0]\n",
       "[0.77716 2.11739 -0.28197; 0.0 1.40951 0.649257; 0.0 0.0 -0.564645]"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Alu = lufact(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Distintas partes de la factorización las puedes accesar con índices especiales"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Array{Float64,2}:\n",
       " 0.0  0.0  1.0\n",
       " 0.0  1.0  0.0\n",
       " 1.0  0.0  0.0"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Alu[:P]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Array{Float64,2}:\n",
       "  1.0       0.0       0.0\n",
       " -0.520989  1.0       0.0\n",
       "  0.203251  0.289228  1.0"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Alu[:L]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Array{Float64,2}:\n",
       " 0.77716  2.11739  -0.28197 \n",
       " 0.0      1.40951   0.649257\n",
       " 0.0      0.0      -0.564645"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Alu[:U]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Podemos calcular la solución de $Ax=b$ del objecto de factorización"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Array{Float64,1}:\n",
       " 1.0\n",
       " 1.0\n",
       " 1.0"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# PA = LU\n",
    "# A = P'LU\n",
    "# P'LUx = b\n",
    "# LUx = Pb\n",
    "# Ux = L\\Pb\n",
    "# x = U\\L\\Pb\n",
    "Alu[:U]\\(Alu[:L]\\(Alu[:P]b))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Más importantemente,* podemos despachar sobre el tipo de `LU` y simplemente resolver el problema por medio de"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Array{Float64,1}:\n",
       " 1.0\n",
       " 1.0\n",
       " 1.0"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Alu\\b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Esto puede ser útil si el mismo lado izquierdo es usado para lados derechos.\n",
    "\n",
    "La factorización también se puede usar para calcular el determinante pues\n",
    "\n",
    "$\\det(A)=\\det(PLU)=\\det(P)\\det(U)=\\pm \\prod u_{ii}$\n",
    "\n",
    "porque $U$ es triangular y el signo está definido por $\\det(P)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.6185193599176367"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "det(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.6185193599176367"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "det(Alu)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### QR\n",
    "Cuando `A` es alta, "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Atall = randn(3, 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Julia calcula la solución de mínimos cuadrados $\\hat{x}$ que minimiza $\\|Ax-b\\|_2$. \n",
    "\n",
    "Esto se puede hacer factorizando\n",
    "\n",
    "```\n",
    "A=QR\n",
    "``` \n",
    "\n",
    "donde $Q$ es unitaria/ortogonal y \n",
    "\n",
    "$R=\\left(\\begin{smallmatrix}R_0\\\\0\\end{smallmatrix}\\right)$ y  $R_0$ es triangular superior.\n",
    "\n",
    "con la factorización QR la norma mínima se puede expresar como \n",
    "\n",
    "\\begin{equation*}\n",
    "\\|Ax-b\\|=\\|QRx-b\\|=\\|Q(Rx-Q'b)\\|=\\|Rx-Q'b\\|=\\left\\|\\begin{pmatrix}R_0x-Q_0'b\\\\Q_1'b\\end{pmatrix}\\right\\|=\\|R_0x-Q_0'b\\|+\\|Q_1'b\\|\n",
    "\\end{equation*}\n",
    "\n",
    "Y entonces el problema se puede reducir a resolver el problema cuadrado $R_0x=Q_0'b$ para $x$.\n",
    "\n",
    "Podemos factorizar QR sobre `Atall` vía"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "ename": "LoadError",
     "evalue": "\u001b[91mUndefVarError: Atall not defined\u001b[39m",
     "output_type": "error",
     "traceback": [
      "\u001b[91mUndefVarError: Atall not defined\u001b[39m",
      "",
      "Stacktrace:",
      " [1] \u001b[1minclude_string\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::String, ::String\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m./loading.jl:522\u001b[22m\u001b[22m"
     ]
    }
   ],
   "source": [
    "Aqr = qrfact(Atall)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Otra característica de la factorización QR es que los tipos `Q` para guardar las matrices unitarias $Q$. Se pueden extraer de tipos `QR` con los índices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Aqr[:Q]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similarmente, la matriz superior triangular $R$ se puede extraer con el índice"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Aqr[:R]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "En este caso se guarda R como una matriz 2x2 en vez de 3x2 porque el último renglón de R está lleno de 0's.<br><br>\n",
    "\n",
    "\n",
    "Aunque la matriz `Aqr[:Q]` se imprime como $3\\times 3$ en el objeto de factorización, en la práctica puede representar la versión delgada también. Así"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Aqr[:Q]*ones(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "funciona y representa a $3 x 2$ matrix por un vector de 2-elementos.\n",
    "\n",
    "Similarmente,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Aqr[:Q]*ones(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "funciona representando la matriz $3x3$ y un vector de 3 elementos.\n",
    "\n",
    "Sin embargo, esto no significa que podemos multiplicar a `Q` por vectores de longitued arbitraria."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Aqr[:Q]*ones(4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La matriz tiene representación interna compacta, entonces indexar sólo hace sentido si uno sabe cómo la factorización guarda datos."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Aqr[:Q][1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "El objeto QRCompactWY `\\` tiene un método para QR y entonces el problema de los mínimos cuadrados es resuelto con"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Aqr\\b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Y si en vez escribimos"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Atall\\b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "En vez de factorizar con QR a `Atall` primero, Julia va a defaultear factorizar *con* pivoteo."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Este default a pivotear la factorización QR le permite a Julia resolver problemas deficientes de rango.\n",
    "\n",
    "Podemos explícitamente escoger pivotear durante la factorización QR (de una matriz singular, por ejemplo) con `Val{true}`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "v = randn(3)\n",
    "# Tomar el producto exterio de un vector consigo mismo nos da una matriz singular\n",
    "singmatriz = v * v'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Aqrp = qrfact(singmatriz,Val{true})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notamos que el tipo que resulta del objecto de Factorization es distinto que antes.\n",
    "\n",
    "`\\` también tiene un método de `QRPivoted` y el problema con rango deficiente es entonces calculado como "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "ename": "LoadError",
     "evalue": "\u001b[91mUndefVarError: Aqrp not defined\u001b[39m",
     "output_type": "error",
     "traceback": [
      "\u001b[91mUndefVarError: Aqrp not defined\u001b[39m",
      "",
      "Stacktrace:",
      " [1] \u001b[1minclude_string\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::String, ::String\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m./loading.jl:522\u001b[22m\u001b[22m"
     ]
    }
   ],
   "source": [
    "Aqrp\\b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Eigendescompisición y los SVDs (Valores de descomposición Singular)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Los resultados de eigendescomposición y de la descomposición singular de valores se guardan en los tipos`Factorization`. Esto también incluye la factorización de Hessenberg y de Schur\n",
    "\n",
    "La eigendescomposición puede ser calculada como"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Asym = A + A'\n",
    "AsymEig = eigfact(Asym)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Los valores y vectores se pueden recoger del tipo Eigen con un índice especial"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "LoadError",
     "evalue": "\u001b[91mUndefVarError: AsymEig not defined\u001b[39m",
     "output_type": "error",
     "traceback": [
      "\u001b[91mUndefVarError: AsymEig not defined\u001b[39m",
      "",
      "Stacktrace:",
      " [1] \u001b[1minclude_string\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::String, ::String\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m./loading.jl:522\u001b[22m\u001b[22m"
     ]
    }
   ],
   "source": [
    "AsymEig[:values]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "ename": "LoadError",
     "evalue": "\u001b[91mUndefVarError: AsymEig not defined\u001b[39m",
     "output_type": "error",
     "traceback": [
      "\u001b[91mUndefVarError: AsymEig not defined\u001b[39m",
      "",
      "Stacktrace:",
      " [1] \u001b[1minclude_string\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::String, ::String\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m./loading.jl:522\u001b[22m\u001b[22m"
     ]
    }
   ],
   "source": [
    "AsymEig[:vectors]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Una vez más, como la descomposición se guarda en un tipo, podemos despatchar sobre esos tipos y explotar un método especializado para cada factorización, e.g. que $A^{-1}=(V\\Lambda V^{-1})^{-1}=V\\Lambda^{-1}V^{-1}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "inv(AsymEig)*Asym"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Julia también tiene una función `eig` que regresa una tupla con los valores y vectores"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "eig(Asym)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "No recomendamos esta versión.\n",
    "\n",
    "La función `svdfact` calcula la descomposición singular de valores"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Asvd = svdfact(A[:,1:2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "y de nuevo `\\` tiene un método para el tipo que permite los mínimos cuadrados por SVD"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Asvd\\b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Existen funciones especiales para proporcionar los valores sólamente: `eigvals` and `svdvals`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Estructuras especiales de matrices\n",
    "\n",
    "La estructura de matrices es muy importante en el álgebra linear. Ésta estructura se le puede hacer explícita a Julia por medio de los tipos compuestos. Ejemplos: `Diagonal`, `Triangular`, `Symmetric`, `Hermitian`, `Tridiagonal` y `SymTridiagonal`. Se han escrito métodos especializados para cada tipo de matriz especial para aprovechar su estructura. Siguen algunos ejemplos:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Creando una matriz Diagonal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Diagonal(diag(A))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Diagonal(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Creando una matriz triangular inferior"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "LowerTriangular(tril(A))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "LowerTriangular(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Creando una matriz simétrica"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Symmetric(Asym)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "SymTridiagonal(diag(Asym),diag(Asym,1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cuando se sabe que una matriz es e.g. triangular o simétrica Julia puede que resuelva el problema más rápido convirtiendo a una matriz especial.\n",
    "\n",
    "Para algunos procedimientos, Julia checa si el input de matriz es triangular o simétrica y lo convierte a tal estructura si lo detecta.\n",
    "\n",
    "Debería notarse que `Symmetric`, `Hermitian` y `Triangular` no copian la matriz original.\n",
    "\n",
    "#### Eigenproblema simétrico\n",
    "La capacidad de Julia para poder detectar si una matriz es simétrica/Hermitian puede influenciar muchísimo sobre qué tan rápido se puede resolver un problema de eigenvalor."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "n = 1000;\n",
    "A = randn(n,n);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Usamos `A` para genera una matriz simétrica `Asym`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Asym = A + A';"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ahora creemos una matriz Asym para simular una matriz simétrica con errores de punto flotante"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Asym_noisy = copy(Asym); Asym_noisy[1,2] += 5eps();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "¿Puede Julia determinar que ambas `Asym` y `Asym_noisy` son matrices simétricas?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "println(\"Is Asym symmetric? \", issymmetric(Asym))\n",
    "println(\"Is Asym_noisy symmetric? \", issymmetric(Asym_noisy))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ahora veamos como el ruido de `Asym_noisy` impacta el tiempo en llevar a cabo una eigendescomposición"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "@time eigvals(Asym);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "@time eigvals(Asym_noisy);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Por suerte, le podeemos proveer información explícita sobre la estructura de la matriz a Julia\n",
    "En este ejemplo, usamos la palabra clave `Symmetric` "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "@time eigvals(Symmetric(Asym_noisy));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Y así nuestros cálculos son mucho más eficientes :)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Un gran problema\n",
    "Usar matrices tridiagonales permite trabajar con problemas potencialmente muy grandes. El siguiente problema no seria posible resolverlo en una laptop si la matriz se tuviera que guardar como tipo `Matrix`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "n = 1_000_000;\n",
    "A = SymTridiagonal(randn(n), randn(n-1));\n",
    "@time eigmax(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Álgebra lineal numérica\n",
    "La manera usual de agregar soporte para álgebra lineal numérica es haciendo un wrapper para subrutinas de BLAS y LAPACK. Para matrices con elementos `Float32`, `Float64`, `Complex{Float32}` ó `Complex{Float64}` esto es lo que hace Julia. Desde hace rato, Julia ha tnido soport para la multiplicación genérica de tipos. Así, cuando uno multiplica matrices de enteros, obtiene una matriz de enteros."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "rand(1:10,3,3)*rand(1:10,3,3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Recientemente, más métodos de álgebra lineal se han añadido a Julia y ahora soporta factorizaciones generales de tipo `LU` y `QR`. Métodos generales para eigenvalores y SVD han sido escritos más recientemente en paquetería externa.\n",
    "\n",
    "En general, la factorización `LU` puede ser calculada cuando los elementos de la matriz se cierraon sobre los operadores `+`, `-`, `*` y `\\`. Por supuesto, la matriz también deben tener rango completo. El método general de factorización `LU` en Julia aplica pivoteo y por lo tanto debe poder soportar `<` y `abs`. Por lo tanto es posible resolver sistemas de ecuaciones de e.g. números racionales como en los ejemplos que siguen."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Para usar números racionales, usamos un doble //:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "1//2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Ejemplo 1: Sistemas racionales lineales de ecuaciones\n",
    "Julia cuenta con números racionales ya instalados. El siguiente ejemplo consta de un sistema lineal de ecuaciones resulto sin promover a tipos de punto flotanto. Puede haber un error de overflow fácilmente al trabajar con racionales, así que usamos `BigInt`s."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Ar = convert(Matrix{Rational{BigInt}}, rand(1:10,3,3))/10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "x = fill(1, (3))\n",
    "b = Ar*x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Ar\\b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "lufact(Ar)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Ejemplo 2: Matriz racional de eigenestructura\n",
    "\n",
    "El siguiente ejemplo muestra como la artimética de matriz racional puede ser usada para calcular una matriz dados los eigenvalores y eigenvectores racionales. Yo he encontrado ésto útil para mostrar ejemplos de sistemas dinámicos lineales."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Array{Rational{Int64},2}:\n",
       " 1//1  -1//2  -1//4\n",
       " 0//1   1//2  -1//4\n",
       " 0//1   0//1   1//4"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "λ1,λ2,λ3 = 1//1,1//2,1//4\n",
    "v1,v2,v3 = [1,0,0],[1,1,0],[1,1,1]\n",
    "V,Λ = [v1 v2 v3], Diagonal([λ1,λ2,λ3])\n",
    "A = V*Λ/V"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercises\n",
    "\n",
    "11.1 ¿Cuáles son los eigenvalores de la Matriz A\n",
    "\n",
    "```\n",
    "A =\n",
    "[\n",
    " 140   97   74  168  131\n",
    "  97  106   89  131   36\n",
    "  74   89  152  144   71\n",
    " 168  131  144   54  142\n",
    " 131   36   71  142   36\n",
    "]\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "11.2 Crea una matriz diagonal de los eigenvalores de A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "11.3 Realiza un factorización de Hessenberg sobre la matriz A. Verifica que `A = QHQ'`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.6.1",
   "language": "julia",
   "name": "julia-0.6"
  },
  "language": "Julia",
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.6.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
